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Abstract. We study the problem of the appropriate choice of the interpolating 
kernel to be used in the evaluation of gradients of functions. Such interpolation 
technique is often used in applications, e.g. it is typical fo Smoothed Particle 
Hydrodynamics (SPH). We propose a minimization procedure for selecting 
kernels in 9ft n , in the class of regular, normalizable, symmetric functions having 
finite moments up to a sufficiently high order; the method is valid when the 
kernel width is position-dependent and allows to recover conservation laws at 
the same order of approximation that SPH, as interpolation technique, has 
when the kernel size is constant. 



1 Introduction 

The aim of this paper is giving a criterion for the choice of kernel functions 
to be used in a particular interpolation problem. The problem is the fitting of 
spatial derivatives of a function / (for example the gradient of /, V/) known 
on a finite set of points in 3? 3 . The results here presented do not depend on 
the way the points are chosen and valid in fully n-dimensional situations. 

The problem of gradient interpolation assumes the form of the right choice 
of interpolating kernels in the so-called Smoothed Particle Hydrodynamics 
(SPH) method (we refer to [1] for an introduction to this numerical method 
and to [2] for some recent improvements). In this scheme the fluid is sampled 
by a set of massive particles which are "fluid representatives": the particle 
positions are considered as moving grid-points for the interpolation of the 
physical quantities characterizing the fluid as well as their spatial derivatives. 
SPH is so a fully Lagrangian method. Analogously to the classic mollifica- 
tion approximation of a sufficiently regular function / with the convolution 
function (/) = /* 0, where * is the convolution operator and is a kernel 
function (whose properties will be specified later), the spatial derivatives of 
/, resumed by V/, are approximated by (V/)i = (V/) * 0. In doing this, 
we take the advantage that, under suitable hypotheses on / and 0, (V/) * 
can be replaced by / * V0. As final step, the integral / * V0 is numerically 
evaluated. Note that the knowledge of / on a set of points suffices to evalu- 
ate (V/)i, and the numerical evaluation of the integral giving / * V0 can be 
done using a Monte-Carlo estimate. The interpolation procedure used in SPH 
follows the scheme we have just sketched, so we shall refer to this scheme as 
"SPH-like" . We start our analysis noting that another approximation for V/ 
can be given by (V/) 2 = V(/*0). As the kernel does not depend explicitely 
on the position, the two approximations (V/)i and (V/)2 are equivalent for 
any regular enough kernel. It is easy to verify that this property is lost when 
the kernel is position-dependent. 

In the applications, to improve spatial resolution the characteristic size of the 
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kernel is actually allowed to depend on the position and the approximation 
(V/)i is implicitely applied. In this case the SPH scheme is no longer mo- 
mentum and energy conservative, [3]. In fact, the use of a position dependent 
kernel introduces, in the approximation (V/)i, a source of error proportional 
to the characteristic size of the kernel while the approximation of the SPH 
scheme is of the second order with respect to it, [1]. 

To overcome this problem, the discrete SPH equations should be modified 
introducing corrective terms as, for example, deeply discussed in [4]. In [4] 
it is shown how the inclusion of the proper V/i terms in the equations of 
motion leads to significant improvement of energy conservation in cases of 
SPH simulations of collisions of polytropic spheres. 

In this paper we approach this problem in another way trying to answer to 
the question: 

• is it possible to reduce unconservativity within the general order of approx- 
imation of SPH without modifying its, appealingly simple, form? 

We believe that the use of the approximation (V/)2 would be desirable for 
this, and we will give, in the Appendix, a strong argument in support. Con- 
sequently, we simply look for kernels such that the exchange between the 
approximation (V/)2 and (V/)i introduces a negligible error. This approach 
has the advantage that no changes on the discrete SPH equations are needed. 

We stress that the kernels proposed and commonly adopted in SPH, have been 
selected either following the spline interpolation theory, [5,6], or minimizing 
the error in the approximation of / with /*</>, [7]. A recent paper, [9], gives 
a quite complete statistical analysis of the quality of kernels of various types 
(i.e. shapes) to reproduce the derivatives of given functions sampled on a one- 
dimensional, even spaced, mesh. Then our analysis is substantially different 
from all the previous ones. In particular we are not interested here in the 
topic of selecting the kernel's width as a function of the particle's position, 
but rather to identify kernels able to keep second order accuracy of SPH as 
interpolation technique even in the case of position-dependent kernel width. 

The plan of the paper is the following: in Section 2 we introduce the kind 
of kernels we deal with; in Section 3 the interpolation procedure is studied 
in detail, and we state a general criterion for selecting kernels; in Section 4 
some kernels are accordingly proposed; in Section 5, a simple numerical test 
is performed, while Sect. 6 is devoted to draw the Conclusions. 
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2 The kernel function 



First of all let us define the kernel to use in the mollification procedure. We 
shall assume the following general form 



where r = (ri,r2, ■ ■■,r n ) is the position vector in 9? n , r = |r|, g(x) is a smooth 
scalar function of the real variable x, with g'(0) = (dg/dx) x=0 = and such 
that the kernel has finite moments up to a sufficiently high order. The nor- 
malization condition / dv — 1 (the integral is extended over all the physical 
space) for any h value is imposed. Notice that these assumptions ensure that 
sup r \f*(f> — f\ oc h 2 , for any sufficiently smooth function /(r). This means that 
the regularized approximating function f * (f) tends uniformly to the "true" 
function / as h goes to zero. The function g defines the shape of the kernel, 
which is particularly important when gradients have to be evaluated. 



3 A criterion for the kernel choice 



In this Section we consider the step of the SPH-like interpolation procedure 
related to the approximation of derivatives of a function. In particular we 
study the difference between (V/) * and V(/ * 0) when depends explicitly 
on the position. As a result we shall derive a first general criterion for the 
right choice of the shape of kernels of the type (1). 

Before starting our analysis we find convenient to specify the notations adopted 
hereafter. We define two differential operators, 



the first one (also referred to as V when no ambiguity occurs) is the usual 
gradient operator in 5J n while the second one is a "partial" gradient operator, 
in the sense that it deals only with the explicit dependence on the spatial 
variables. 

As we have sketched in the Introduction, the basic idea behind the SPH- 
like method is the approximation of the gradient of a quantity basing on the 



0(r, h) oc h n g(r/h) 



(1) 
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knowledge of / on the particle positions, only. This idea relies on that 



(V/)*0 = /*V0 = V(/*0), 

these identities being true only under appropriate boundary conditions for / 
and and if h = const. When h = h(r) then integration by parts leads to 
V(/ * 0) = / * <90 and 

(V/) * = (V/)! = J V r ,/(r')0(r - r ', h(r))dr> = 
= V r J /(r')0(r - r', h(r))dv' - 5 = V(/ * 0) - 5 = (V/) 2 - 5, (3) 



where 5 depends on / and h. In this case, the normalization condition on the 
kernel implies 

J <90(r - r', h(r))dr' = 0, J J^0(r - r', h(v))dr' = 0. (4) 



Our aim is now to find how making 5 as small as possible independently of /. 
Taking into account (1), let us evaluate 5: 

5 = J /(r')VMr)^0(r-r',/i(r))dr'oc 



In order to recover the equality (V/) * = V(/ * 0), independently of /, the 
function g in (5) should satisfy, in K + , the ordinary differential equation 

ng{x) + xg\x) = 0. (6) 



The non-trivial solution of this ODE is given by g oc x~ n , which leads to a 
non-normalizable kernel. This means that, with a choice of g in the set of 
normalizable kernels, 5 cannot vanish for all functions /. 

Thus, a possibility is to minimize 5 via its expansion in terms of the parameter 
h, looking for conditions that allow to eliminate at least some low order terms 
on h. 

Performing a Taylor expansion of the function / up to the second order around 
the point r, and taking into account the symmetry of the kernel and the 
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normalizing conditions, it is easy to verify that 

6 = ^VMr)A/(r) J (n - A) 2 J>(r - r', h(v))dr' + ^)0(h 2 ), (7) 

where ip is a function of the position and A is the Laplace operator. Now, we 
look for kernels that minimize the absolute value of the integral in (7). Let 
us work out the case n > 3 (the simpler cases n = 1 and n = 2 must be 
handled in a slight different way). Making the substitution z = (r — r')/h(r) 
the integral in (7) can be written as 

h(r)jzi(ng(\z\)+g>(\z\)\z\)dz. (8) 

To perform the integral in (8) we introduce the set of cylindrical variables 
p, #3, ...,#„), thus getting, after simple manipulations of variables: 

^ oo oo 

-(n- l)w n _x J dp J \p\ n - 2 zf 

— oo — oo 

where u n is the volume of the n-dimensional unitary sphere. Moreover, using 
polar coordinates (a, (5) on the plane (z 1 ,p), the expression (9) becomes 

2-7T OO 

-{n - l)w n _i y sm 2 (/3)|cos(/3)|rf/3 ^ a n+1 (n#(a) + g\a)a)da; (10) 
o o 

the normalization condition, now, is 

_^ 2-7T OO 

-(n-l)u n ^J \sin((3)\d(3 J g(a)a n ~ 1 da= 1. (11) 
o o 

It is clear that the approximation of (V/) * with V(/ * (f>) improves of, at 
least, one order with respect to h if the normalized kernel has a shape such 
that the integral over a in (10) is zero. This improvement makes the SPH-like 
interpolation with variable kernel size of the same order (the second) of the 
classic SPH interpolation with constant h (see Appendix) . 

Thus, we get a criterion to select kernels: 

kernels of type (1) should be chosen selecting g in the class of differentiable 



ng 



2 + P 2 



2 + P 2 



dzi, 
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functions such that 

oo oo 

J a n+1 {ng{a) + g'(a)a) da = , 0< J a^g (a) da < oo. (12) 
o o 

Remark 1: it is easily shown that also the cases n — 1 and n = 2 lead to the 
same constrain (12) for the selection of the kernel. Remark 2: the stated crite- 
rion is of course useless if V/i(r) is imposed to be zero. Remark 3: this analysis 
may be iterated by evaluating higher order terms in the Taylor expansion of 
the interpolated function /, and imposing them to be zero. 

As final, relevant, consideration we notice that the shape of our proposed 
'optimal ' kernel depends on the dimension, n, of the space ; of course this 
should be taken into account in the practical applications. 



4 Some examples and applications 



We want, now, to apply the stated criterion to select kernels (in some specific 
class of functions) in the most common three-dimensional case. The most 
used kernels in SPH numerical applications derive from the theory of spline 
functions, [5,6]. Thus, we consider the class of functions (which corresponds 
to the third-order spline functions class) defined as: 



g(r) = 



ar 3 + br 2 + cr + d 
a'r 3 + b'r 2 + c'r + d' 




< r < 1 

1 < r < 2 
r > 2, 



(13) 



continuous up to the second derivative and such that: 

(1) g(0) = 1 and g(2) =0; 

(2) the first derivative vanishes on the boundary of the set [0, 2]. 

Following the criterion discussed in Section 3, the selected normalized 3-D 
kernel in this class is found to be 



<j>ir,h) = 



15 



11527r/i 3 



f (r/h) 2 (171 (r/h) - 321) + 172 
((r/h) - 2) 2 (107 - 85 (r/h)) 




< r/h < 1 

1 < r/h < 2 (14) 
r/h > 2. 
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Fig. 1. The kernel (14) selected with our method (solid line) is compared with the 
kernel (17) (dashed line). 

Another kernel often used in the applications belongs to the super-Gaussian 
class of kernels: 

g (r) = (ar 2 + 6)e~ r2 , r > 0. (15) 



When applied to this class of functions, our criterion yields (in 3-D) 

^h) = J^(l-(r/hf)e-^\ (16) 



We note that the kernel (14) differs from the kernel belonging to the same 
class and suggested by Monaghan and Lattanzio [6] (see also Figure 1), i.e.: 



</>(r,h) 



irh 3 



l-\{r/hf + \{r/hf 
\ (2 - {r/h)f 




< r/h < 1 

1 < r/h < 2 (17) 
r/h > 2, 



On the other hand the kernel (16) is exactly the same super-Gaussian kernel 
suggested by Gingold and Monaghan [7]. 

Finally, note that a kernel which is somewhere negative is compatible with our 
analysis, for we never made assumptions on the sign of g in (1). Nevertheless, 
a non-negative kernel should be desirable whether the interpolated function is 
thought to have some specific physical meaning (for example when it represents 
a mass density). A possible way to overcome this problem could be the use of 
two different kernels to represent the physical quantities and their gradients 
separately, but this point should deserve a deeper analysis. 
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5 A numerical test 



The problem of having a good interpolation of derivatives occurs often in the 
applications and, even if it does not guarantee the quality of SPH evolutive 
simulations, surely helps in fluid-dynamical numerical simulations. 

Of course, when violent shock fronts appear in fluid-dynamical simulations, to 
keep energy and momentum conservations at a high level of precision requires 
both a large number of particles in around the front (i.e. high resolution) and 
modification of SPH equations by mean of the inclusion of the relevant V/i 
term, as shown in [4]. Actually, we could check this in an SPH simulation 
of the classic case of 1-D shock tube calculation, after selecting the best 1-D 
kernel in the class of cubic splines, following our, previously described, recipe. 
The treatment of shocks improves with our kernel, but the most significant 
improvement come actually from higher resolutions and V/i terms inclusion. 

We shall use here our selected kernel (14) in a numerical example of function 
gradient interpolation, comparing our results with the analogous ones obtained 
using the kernel (17) suggested by Monaghan and Lattanzio [6] (of course, the 
3-D case is by far the most important in practical applications). To this aim, 
we evaluate numerically the gradient of a function /, supposed known on a set 
of positions, by SPH-like interpolation. We will perform the test using kernels 
whose widths strongly depend on the position. To this end we will set the 
points on a non-uniform grid, and will fix the kernel widths so that only a 
fixed number of points will be "touched" by any kernel. 

The integral (V/) * <fi ~ V(/ * 0) is then approximated, on the positions r^, 
by the Monte-Carlo estimate 



where N is the number of positions rj over which / is known and p is the 
number density of the points Tj. As / is assumed to be density of the points, 
the Monte-Carlo estimate reduces to the simpler form 




i = 



1 N 

V r /W~^£ (9 r 0(r-r„Mr))) r=ri 



(18) 



In this case, the function / itself is approximated as follows 



1 



N 



f(n) 



N 



%*))• 



(19) 
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For the sake of comparison we perform the numerical evaluations using both 
the kernels (14) and (17). 

As "test" function, we choose the function f(x,y,z) = (2|x — 1/2|)~ 1 / 2 in the 
cube < x, y,z < 1, assuming / as point density, too. Being / diverging as 
x — > 1/2, a lot of points r,- will fall in the neighbourhoud of x — 1/2. The 
exact gradient is V f(x,y, z) = ((x — l/2)/(2y / 2|a: — 1/2 1 5 ), 0, 0), singular in 
x = 1/2. Taking into account the meaning of / as a density, the points are 
displaced on a regular grid. It is easy to verify that the i x — th position is 
given by the root of the expression 1 + \pl • (x — 5)/\J\x — 0.5 1 — 2 ■ i x /n x , 
where 1 < i x < n x and n x is the number of grid points on the cube edge 
in the x direction. The setting on the other coordinates is trivial. The total 
number of grid points is so n x ■ n y ■ n z . In the test presented here n x = 60 and 
n y = n z = 50. The fit to the function and to its gradient have been obtained 
applying formulas (19) and (18) respectively. 

The kernel sizes h are chosen such that for each point the sums (19) and 
(18) have only 2.5% of non-zero terms, then the size of the kernel depends 
strongly on the position. That is, the width h changes very rapidly in space. 
This is surely a situation for which 5, in equation (3), differs from zero: this 
fact justifies the choice of this particular /. 

In Figure 2 the fits obtained with the two cited kernels, along the line y = z = 
1/2, are presented. Our kernel seems to fits quite better either the function / 
either the x derivative. 



Fig. 2. In Panel (a) the interpolation of the funcion / (solid line) obtained using our 
kernel (dotted line) and Monaghan and Lattazio's (dashed line) are shown; Panel 
(b) refers to df/dx. 

As a final consideration, we note that, from a theoretical point af view, the 
results of this test should depend mainly on the specific dependence of the 
kernel's widths, h = h(r), with respect to the position, rather than on the 
particular function, /, chosen for the test. 
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6 Conclusions 



In this paper we have studied an interpolation problem of great interest in 
practical applications. The problem is the approximation of the derivative of 
a function, known on a set of points, by a kernel interpolating method. Given 
a suitable kernel function 0, the idea of the method consists in approximating 
the derivatives of a function /, and so its gradient, with its mollified version 
(V/)i = (V/) *(/) = /* V0, and then evaluating numerically the last convolu- 
tion integral. This procedure fails as the kernel size depends on the position, 
i.e. as the operation "*" cannot be interpreted as a convolution. In this case, 
we stated a rule to choose the kernel. The kernel is chosen such to control 
the difference between (V/)i and the other possible mollification defined as 
(V/)2 = V(/ * 0) up to a given order in h. This simple criterion is indepen- 
dent of the numerical method applied in the integration. We have given some 
examples of kernels selected applying our criterion and discussed results in a 
typical case. 



Appendix 



Let us consider the evolutive equations for a compressible fluid written in 
conservative form, [8]: 



Pt 

Qit + V r 
e t + V r 



V r • q = 



1,2,3; 



(20) 



0, 



where q = pv is the density of momentum, P is the pressure, e/p — v 2 /2 + e 
is the total energy per unit mass and e is the internal one. 

The usual way to get the discrete SPH equations consists in the convolution of 
each equation with a kernel and then in its approximation through Monte- 
Carlo integrations, [1]. We shall consider here the first step only. We begin 
considering a kernel of the type (1) where h is a constant. The set of equations 
(20) becomes: 



(<li)t 



(P)t + V r • (q) = 
V r -((qf)) + V ri (P) 



1,2,3; 



(21) 



(e) t + V r 



([(e + P)* 



0, 



where (p) = J p(r')0(r — r', h)dr', and so on. Let us study the conservation of 
the mollified quantities. It is straightforward to verify that ^ f(p)dr — and 
f(p)dr — J p dr, and similar for the other quantities. The conservation of the 
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mass derives directly by the equation (V • q) * <fi = V • (q * <$) which means 
(see the Introduction) (V • q)i = (V • q) 2 . 

When h = h(r), then | J(p)dr = 0(h) and so on (see Section 3); that is the 
system is not longer conservative and the error derives from the exchange of 
(V • q) * <p with V • (q* 0). In other words, the error derives from the exchange 
of (V-q)i with (V-q)2- The conclusion is that when the approximation (V/)2 
cannot be applied, the conservations are lost. 

The use of kernels selected on the basis of the criterion stated in Section 3 
ensures -| f{p)dr = 0(h 2 ), etc., which is the order of approximation of the 
SPH method, [3]. 



References 

[1] J.J. Monaghan, Smoothed Particle Hydrodynamics, Ann. Rev. of Astron. and 
Astroph. 30 (1992) 543-574. 

[2] P.W. Randies and L.D. Libersky, Smoothed particle hydrodynamics: some 
recent improvements and applications, Comput. Methods Appl. Mech. Engrg. 
139 (1996) 375-408. 

[3] G.V. Bicknell, The equations of the motion of particle in Smoothed Particle 
Hydrodynamics, SIAM J. Sci. Stat. Comput., 12-5 (1991) 1198-1206. 

[4] R.P. Nelson and J.C.B. Papaloizou, Variable smoothing lengths and energy 
conservation in smoothed particle hydrodynamics, Mon. Not. of the Roy. Astr. 
Soc, 270 (1994) 1. 

[5] J.J. Monaghan, Extrapolating B -splines for interpolation, J. Comp. Phys. 60-2 
(1985) 253-262. 

[6] J.J. Monaghan and J. Lattanzio, A refined particle method for astrophysical 
problems, Astron. and Astroph. 149 (1985) 135-143. 

[7] R.A. Gingold and J.J. Monaghan, Kernel estimates as a basis for particle 
methods, J. Comp. Phys. 46 (1982) 429-453. 

[8] L.D. Landau and E.M. Lifshitz, Fluid Mechanics (Pergamon, 1959). 

[9] D. A. Fulk and D. W. Quinn, An analysis of 1-D smoothed particle 
hydrodynamics kernels, J. Comp. Phys. 126 (1996) 165-180. 



12 



